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1. Motivation and objectives 

Free shear flows, like those of mixing layers, are encountered in aerodynamics, in 
the atmosphere, and in the ocean as well as in many industrial applications such 
as flow reactors or combustion chambers. It is, therefore, crucial to understand the 
mechanisms governing the process of transition to turbulence in order to predict 
and control the evolution of the flow. Delaying transition to turbulence as far 
downstream as possible allows a gain in energy expenditure while accelerating the 
transition can be of interest in processes where high mixing is desired. Various 
methods, including the use of polymer additives, can be effective in controlling fluid 
flows. 

The drag reduction obtained by the addition of small amounts of high polymers 
has been an active area of research for the last three decades. It is now widely 
believed that polymer additives can affect the stability of a large variety of flows and 
that dilute solutions of these polymers have been shown to produce drag reductions 
of over 80% in internal flows and over 60% in external flows under a wide range of 
conditions. (Berman 1978, Sellin 1985 and Sellin & Moses 1989) 

The major thrust of this work is to study the effects of polymer additives on 
the stability of the incompressible mixing layer through large scale numerical sim- 
ulations. In particular, we focus on the two-dimensional flow and examine how 
the presence of viscoelasticity may affect the typical structures of the flow, namely 
roll-up and pairing of vortices. 

2. Accomplishments 

2.1 Problem definition 

The flow is examined in a reference frame moving with the average velocity. In 
such frame, the flow is characterized by the upper free-stream velocity u Q and the 
momentum thickness of the mixing layer 6. We used a vorticity-streamfunction 
formulation for Cauchy’s momentum equation, 

Dv 

P— = -Vp+V-T (1) 

This equation is closed through evolution equations relating the stress tensor to the 
shear rate tensor. In all the subsequent analysis, the stress tensor is written as the 
sum of two terms (Larson (1988) and Bird et al. (1987)): 

T = T* + T p = r/ s 7 + ri p a = t)[k7 + (1 - *)a] (2) 
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The first term corresponds to the contribution of the Newtonian solvent and is 
proportional to the shear rate tensor with t} s the solvent viscosity. The second 
term represents the polymeric contribution to the stress, and is proportional to the 
tensor a with rj p the polymeric contribution to the shear viscosity. The parameter 

jo yi 

k = — = — may vary between 0 and 1. When k — 1, the field equations 

V* + Vp V 

and the constitutive equations can be decoupled, and the problem reduces to that 
of purely Newtonian flow. 

The tensor a satisfies the appropriate rheological equation that can be of differ- 
ential or integral form (Bird 1967). In the present study, we used two rheological 
models, the Oldroyd-B model and the FENE-P model. 

2.2 Rheological models 

In the Oldroyd-B model, the polymer stress a satisfies the upper convected 


Maxwell equation: 


6a 

A- + a = T 

(3) 

where: 

II 

* Ics 

da. _ . _ 

— — \- v ■ Va — Vt> * a — a • Vv 
ot 

(4) 


is the upper- convected derivative of a, and A is the relaxation time of the poly- 
mer. This model gives a reasonably good qualitative description of dilute polymer 
solutions, but unfortunately, it gives rise to a steady state elongational viscosity 
that diverges at a finite elongational rate. This unlikely behavior results from the 
infinite extensibility of the linear Hookean spring used to model the polymer. In 
order to avoid this problem, a nonlinear spring based on Warner law can be used 
to describe the finite extensibility of the polymer, leading to the FENE-P model. 

U/pD\ 

This model is best formulated in terms of the tensor B = — — — where (RR) is 

the configuration tensor, H the spring constant, k the Boltzmann constant, and T 
the absolute temperature. The tensor B satisfies the following equation: 

(5) 


ZB + A^ = 1 
ot 


J? 2 /rB 

In Eq. (5), I is the unit tensor and Z = (1 — ( — ^)) -1 = (1 — ) -1 . The parameter 

Rq V 

HR 2 

b = — ■ , where R 0 represents the maximal possible extension of the polymer, is 

a measure of the extensibility of the polymer chain. An equivalent formulation of 
Eq. (5) in terms of the tensor a can be obtained using the transformation: 


a = 


ZB -I 
A 


(6) 
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2.3 Scaling and parameters 

Using u Q and 6 as the reference velocity and the reference length, respectively, 

p6u 0 6u 0 

the flow is characterized by the Reynolds number, Re = = , where v is 

T] u 

Au 0 

the kinematic viscosity of the fluid, and the Weissenberg number, We = is a 

dimensionless measure of the relaxation time of the polymer. The elasticity number 

E = ^ i s often used to characterize the elasticity of the fluid. In addition 

Re 

to Re and We, k and, in the case of the FENE-P model, b are the other model 
parameters. 

2.4 Numerical method 

The simulations reported in this study were performed by solving the vorticity 
equation: 

. d k —o-i r d d , (1 - /c) r/ d 2 d 2 . d 2 

-r — V 2 ]iJ — “ U— + 17 — [(*o ^ ~2k~~2 )^12 q * 

l dt Re 1 dx dy Re dx 2 dy 2 dxdy 


(<322 — a l 1 )] (7) 


coupled with the appropriate stress equations. In the present study, we are inter- 
ested in the forced, temporally growing mixing layer. The initial flow is composed 
of the viscously spreading tarth vorticity profile and the corresponding base-state 
polymer stress, seeded with the wave that, according to linear stability analysis 
(Azaiez &; Homsy 1993), has the largest growth rate. 

The dynamical equations are solved using a pseudo-spectral method in which the 
flow variables are expanded in a modified Hartley series (Zimmerman & Homsy 
(1991)). The resulting set of ordinary differential equations is advanced in time 
using an operator splitting algorithm (see e.g. Tan & Homsy (1987)). In addition 
to those with the spectral code, a few simulations have been conducted using a 
finite difference scheme second order accurate in space and in time. A description 
of the scheme can be found in the paper by Orlandi et al. (1992). The results 
obtained from these two codes were always in total concordance. The codes were 
validated by comparing with the linear stability results (Azaiez & Homsy (1993)) 
and by checking that they reproduced the same results as the Newtonian code when 
we set k = 1. 

A typical run for the roll-up of the non-Newtonian fluid required 128x128 spectral 
modes and a time step At = 0.04. This resolution gives satisfactory results at 
moderate values of the Weissenberg number and was refined for large values of We. 
Throughout this study, the value of the parameter k is fixed to 0.5. 


3. Results 

3.1 The Oldroyd-B model 

The vorticity and stress equations for the Oldroyd-B model have been solved nu- 
merically for various values of W e, and for Re =100. For small values of We ( We ~ 1 ), 
the flow does not show any noticeable changes from the Newtonian case. Numer- 
ical simulations at moderately high We (We ~ 10) developed an instability that 
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T = 30.0 



THE NEWTONIAN FLUID. Re=100. 





THE FENE-P MODEL. Re=100, We=50, b=5. 


FIGURE 1. Vorticity contours for Re=100 
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T = 20.0 T = 30.0 

Minimum Contour= 5.000E-01 Minimum Contours 5.000E-01 

Maximum Contour= 1.310E+00 Maximum Contour= 1.260E+00 



T = 40.0 T = 50.0 

Minimum Contours -5.000E-01 Minimum Contour= -3.000E-01 

Maximum Contour= 1.100E+00 Maximum Contours 1.000E+00 



T = 60.0 T = 70.0 

Minimum Contours 4.000E-01 Minimum Contour 4.000E-01 

Maximum Contour= 1.000E+00 Maximum Contours 1.000E+00 



T = 80.0 T = 90.0 

Minimum Contours 4.000E-01 Minimum Contours 4.000E-01 

Maximum Contours 1.000E+00 Maximum Contours 1.000E+00 


Figure 2. (Bll -B22) contours for FENE-P model. Re=50, E=0.5, and b=5. 
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lead to the divergence of the code. We examined the origin of this instability by 
looking at the evolution of the different terms in the polymer stress equations. The 
analysis of these terms showed that the instability is associated with a deficiency 
of the Oldroyd-B model that allows stresses to grow indefinitely. The instability 
starts to develop first in the braid regions where the product of the Weissenberg 
number and the dimensionless local extensional rate exceeds unity. In these regions, 
and due to high extensional rates, the chain is stretched rapidly, and because of its 
large relaxation time associated with the high We, it is prevented from coiling up 
as quickly as its stretching. As a consequence, the chain gets extended indefinitely 
and the stresses grow exponentially. The intense build-up of the stresses ultimately 
leads to the divergence of the numerical code. 

3.2 The FENE-P model 

Unlike the Oldroyd-B model, the FENE-P model does not allow infinite extension 
of the spring used to model the polymer, and as we have seen, the maximal extension 
of the spring is characterized by the parameter b. The viscoelastic mixing layer has 
been successfully simulated for various values of the three parameters, Re,E , and 
b . In what follows, we describe results for the two mechanisms of instability of the 
two-dimensional mixing layer, namely the roll-up and the pairing of the flow. 

3.2.1 Roll-up 

We explored values of the Reynolds number between 50 and 400, varied the 
Weissenberg number between 20 and 200, and examined values of b between 1 and 
20. Figure 1 shows a time sequence of the roll-up of the flow for the Newtonian case 
and for the FENE-P model. As it has been experimentally documented (Riediger 
1988), we observed a trend for smaller values of the minimal (negative) vorticity in 
the case of the viscoelastic flow, as well as a tendency for the vortex structures to 
be more compact and to have longer life times than in the Newtonian fluid. The 
global structure of the flow as well as the roll-up time are basically the same for both 
fluids. However, the local distribution of the vorticity is affected by the presence of 
the viscoelasticity in the flow with high gradients tending to appear in some parts of 
the flow, namely in the braids. The evolution of the absolute value of the minimal 
vorticity at various streamwise locations confirms the conclusion of the tendency to 
have more spanwise vorticity remaining in the braid region. 

The examination of contours of the first normal polymer stress (Bn — B 22 ) showed 
that there is a spatial correlation between the regions of intensification of the vor- 
ticity gradients and those where there is a build-up of the first normal polymer 
stress. 

In order to understand the reasons why the global structure of the roll-up remains 
unchanged, we examined the evolution of the polymer stresses in connection with 
that of the vorticity (Figure 2). This study revealed that the first normal stresses 
reach a quasi- steady state characterized by the absence of any extensional forces 
and a balance between shearing forces and the polymer relaxation stresses, and it 
is interesting to note the spatial relationships between vorticity and normal stresses 
that characterize this quasi-steady state. After the first stage of roll-up, most of the 




T = 40.0 T = 40.0 

Minimum Contours -3.300E-01 Minimum Contour -3.900E-01 

Maximum Contour= -7.000E-02 Maximum Contours -7.000E-02 



T = 60.0 T = 60.0 

Minimum Contour= -2.900E-01 Minimum Contours -3.400E-01 

Maximum Contours -7.000E-02 Maximum Contours -7.000E-02 



T = 80.0 T = 80.0 

Minimum Contours -2.600E-01 Minimum Contours -3.000E-01 

Maximum Contours -7.000E-02 Maximum Contours -7.000E-02 



T =100.0 T =100.0 

Minimum Contours -2.400E-01 Minimum Contours -2.700E-01 

Maximum Contours -7.000E-02 Maximum Contours -7.000E-02 



T =120.0 T =120.0 

Minimum Contours -2.200E-01 Minimum Contours -2.500E-01 

Maximum Contours -7.000E-02 Maximum Contours -7.000E-02 


NEWTONIAN FLUID, Re=50 FENE-P MODEL Re=50, E=1 , b=5 

Figure 3. Vorticity contours for FENE-P model. Re=50, E=l, and b=2. 
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vorticity is entrained into the vortex core with only little remaining in the braids, 
while the stress fields are highest in regions of highest shear, which tend to lie in 
nearly irrotational regions surrounding the vortex core. This feature may be key to 
explaining why the global flow structure is unaffected by viscoelasticity. 

Simulations at various values of the parameters E and 6 show that the estab- 
lishment of the steady state for (Bn - B 22 ) is, at least for the range of E and b 
examined, a common trend of the evolution of the polymer stresses in free shear 
flows and is insensitive to the value of E or b. A detailed discussion of this steady 
state is found in the article by Azaiez & Homsy (1994) 

8.2.2 Pairing 

We examined the pairing of the flow for Re = 50 and various values of E and b. In 
the bulk of the simulations, we used 256 grid points in the streamwise direction and 
128 grid points in the transverse direction and a time step At = 0.02. Simulations at 
large E (E ~ 2) required a finer resolution and we used 256x256 spectral modes and 
a time step At = 0.01. Figure 3 shows a time sequence of the vorticity contours for 
the Newtonian fluid with Re=50 and the non-Newtonian fluid with Re = 50, E = 1 
and 6 = 5. In the early stage, the flow shows the same trends for intensification 
of the vorticity that we have encountered in the case of the single roll-up. Later 
the two vortices start their orbital motion, with the tendency for a slightly faster 
rotational motion in the case of the viscoelastic fluid. Note that during roll-up as 
well as pairing, the vortices are more diffuse in the case of the Newtonian fluid and 
that the maximal absolute value of the vorticity over the whole flow is larger for 
the viscoelastic fluid. We attribute the faster rotation of the two parent vortices 
around each to the vorticity gradients that develop in the braids during the roll-up 
of the flow. These gradients lead to a stronger outer field between the two parent 
vortices that enhances the mutually induced rotational motion of the two vortices 
(Azaiez & Homsy (1994)) 

4. Conclusion 

In the present study, we examined the instability of the plane, incompressible, 
non-Newtonian mixing layer, focusing on simulations with high Re and W e. Numer- 
ical simulations using the Oldroyd-B model developed an instability for moderate 
We. We examined the origin of this instability by looking at the evolution of the dif- 
ferent terms in the stress equations, which showed that the instability is associated 
with a deficiency of the Oldroyd-B model that allows stresses to grow indefinitely. 
The instability starts to develop first in the braid regions where the product of 
the Weissenberg number and the local extensional rate is larger than one is larger 
than one. The unbounded growth and intense build-up of the stresses ultimately 
leads to the divergence of the numerical code. Most of the numerical simulations 
have been performed with the FENE-P equations which revealed to be the most 
appropriate model for the simulation of free shear flows at high elasticity. These 
simulations showed that, for the range of parameters examined, the global structure 
of the flow as well as the roll-up and pairing times are unchanged from the New- 
tonian case. However, local vorticity intensifications associated with the build-up 
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of normal stresses have been observed in the braids as well as in the vortex core. 
As it has been experimentally documented (Riediger 1988), we observed a trend for 
smaller values of the minimal (negative) vorticity in the case of the viscoelastic flow 
as well as the tendency for the vortex structures to be more compact and to have 
longer life times than in the Newtonian fluid. 

The examination of the evolution of the first normal stresses revealed a very 
interesting steady state characterized by the absence of any extensional forces and 
a balance between shearing forces and the polymer relaxation stresses. 
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